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point set in 1R 3 . Each shape is a well-defined polytope, derived from the Delaunay triangulation 
of the point set, with a parameter a £ 1R controlling the desired level of detail. An algorithm is 
presented that constructs the entire family of shapes for a given set of size n in time 0(n 2 ), worst 
case. A robust implementation of the algorithm is discussed and several applications in the area 
of scientific computing are mentioned. 

Categories and Subject Descriptors. F.2.2 [Analysis of Algorithms and Problem Complex- 
ity]: Nonnumerical algorithms and problems — geometrical problems and computations; G.4 [Mathe- 
matics of Computing]: Mathematical Software — reliability and robustness; 1.2.10 [Artificial Intel- 
ligence]: Vision and Scene Understanding — representations, data structures, and transforms; shape; 
1.3.4 [Computer Graphics]: Graphics utilities — application packages; graphics packages; 1.3.5 [Com- 
puter Graphics]: Computational Geometry and Object Modeling — curve, surface, solid, and object 
representations; geometric algorithms, languages, and systems; J. 2 [Computer Applications]: Physical 
sciences and engineering. 

General Terms. Algorithms, software, reliability. 

Additional Key Words and Phrases. Computational graphics, geometric algorithms, scientific vi- 
sualization, scientific computing; three-dimensional space, point sets; polytopes, simplicial complexes, 
Delaunay triangulations; robust implementation, simulated perturbation. 

ACM Transactions on Graphics, 13(l):43-72, 1994. 

1 Research of both authors is supported by the National Science Foundation under grant CCR-8921421. Work 
by the first author is also supported by the Alan T Waterman award, grant CCR-9118874. Any opinions, findings, 
conclusions, or recommendations expressed in this publication are those of the authors and do not necessarily reflect 
the view of the National Science Foundation. 

2 Address: Department of Computer Science, University of Illinois at Urbana-Champaign, 1304 W Springfield 
Ave, Urbana, IL 61801, USA. 



1 Introduction 



The geometric notion of "shape" has no associated formal meaning. This is in sharp contrast 
to other geometric notions, such as diameter, volume, convex hull, etc. The goal of this paper 
is to offer a concrete and formal definition of shape that can be computed and applied. It is 
not supposed to possibly cover the entire range of meanings the term "shape" carries in our 
contemporary language, even if restricted to geometric contexts. Nevertheless, it is sufficiently 
flexible to facilitate a wide range of applications including the study of molecular structures and 
the distribution of galaxies in our universe (see section 7). 

More specifically, the topic of this paper is the definition and computation of the shape of a finite 
point set in three-dimensional Euclidean space, 1R 3 . Intuitively, we think of the set as a cloud of 
points and we talk about the shape of this cloud. A peculiar aspect of the common usage of the 
word "shape" is that its meaning varies with the degree of detail intended. This aspect will be 
reflected by defining a one-parametric family of shapes ranging from "fine" and "local" to "crude" 
and "global". 

A fair amount of related work has been done for point sets in 1R 2 } and some for point sets in 1R 3 . 
Jarvis [23] was one of the first to consider the problem of computing the shape as a generalization 
of the convex hull of a planar point set. A mathematically rigorous definition of shape was later 
introduced by Edelsbrunner, Kirkpatrick, and Seidel [14]. Their notion of a-shapes is the two- 
dimensional analogue of the spatial notion described in this paper. Two-dimensional a-shapes are 
related to the dot patterns of Fairfield [17, 18] and the circle diagrams used in bivariate cluster 
analysis (see for example Moss [33]). Different graph structures that serve similar purposes are 
the Gabriel graph [31], the relative neighborhood graph [39], and their parameterized version, the 
/3-skeleton [26]. 

For 1R 3 } Boissonnat [1] suggested the use of Delaunay triangulations in connection with heuristics 
to "sculpture" a single connected shape of a point set. Our concept of shape is more general and 
mathematically well defined. More recently, Veltkamp [40] also generalized the above-mentioned 
two-dimensional graph structures to three-dimensions, calling them 7-graphs. Finally, note the 
superficial similarity between a-shapes and isosurfaces in 1R 3 . The latter is a popular concept in 
volume visualization (see for example [8,29]). 

The outline of this paper is as follows. A formal definition of a-shapes, along with illustrations, is 
presented in section 2. Geometric concepts related to a-shapes are discussed in section 3. These are 
a- hulls, a-diagrams (also known as space-filling diagrams), a-complexes, Delaunay triangulations, 
and Voronoi diagrams. A combinatorial analysis of a-shapes is offered in section 4. Using Delaunay 
triangulations, it is fairly easy to compute a-shapes in 1R 3 . The resulting algorithm is sketched in 
section 5. Given a set of n points in 1R 3 } it constructs a convenient implicit representation of the 
family of all a-shapes in time 0(n 2 ) } worst case. This algorithm has been implemented by the 
second author of this paper. In section 6 we report on essential aspects of the implementation, 
such as its data structures, how it achieves robustness, and how it performs in practice. Section 7 
discusses some application problems that might benefit from the use of a-shapes. Finally, section 8 
considers possible extensions of the material presented in this paper. 



1 



2 Alpha Shapes in Space 



This section gives an intuitive description as well as a formal definition of three-dimensional a- 
shapes. Both are supported by illustrations that show point sets with sample members of their 
a-shape family. The beauty and elegance of the concept of an a-shape will hopefully be obvious 
after reading section 3 where relationships to other natural geometric concepts are revealed. 

Intuitive Description. Conceptually, a-shapes are a generalization of the convex hull of a point 
set. Let S be a finite set in 1R 3 and a a real number with < a < oo. The a-shape of S is a 
polytope that is neither necessarily convex nor necessarily connected. For a = oo, the a-shape 
is identical to the convex hull of S. However, as a decreases, the a-shape shrinks by gradually 
developing cavities. These cavities may join to form tunnels, and even holes may appear (see 
figure 1). 

Intuitively, a piece of the polytope disappears when a becomes small enough so that a sphere 
with radius a, or several such spheres, can occupy its space without enclosing any of the points 
of S. Think of 1R 3 filled with styrofoam and the points of S made of more solid material, such 
as rock. Now imagine a spherical eraser with radius a. It is omnipresent in the sense that it 
carves out styrofoam at all positions where it does not enclose any of the sprinkled rocks, that is, 
points of S. The resulting object will be called the a-hull (see section 3). To make things more 
feasible we straighten the surface of the object by substituting straight edges for the circular ones 
and triangles for the spherical caps. The obtained object is the a-shape of S (see figure 2). It 
is a polytope in a fairly general sense: it can be concave and even disconnected, it can contain 
two-dimensional patches of triangles and one-dimensional strings of edges, and its components can 
be as small as single points. The parameter a controls the maximum "curvature" of any cavity of 
the polytope. 

General Position. Throughout this paper we assume that the points of S are in general position. 
For the time being, this means that no 4 points lie on a common plane; no 5 points lie on a common 
sphere; and for any fixed a, the smallest sphere through any 2, 3, or 4 points of S has a radius 
different from a. The general-position assumption will later be extended when convenient (see 
section 5.3). It simplifies forthcoming definitions, discussions, and algorithms, and is justified by a 
programming technique known as SoS [15]. This method simulates an infinitesimal perturbation 
of the points on the level of geometric predicates and relieves the programmer from the otherwise 
necessary case analysis (see section 6.2). 

Formal Definition. For < a < oo, let an a-ball be an open ball with radius a. For complete- 
ness, a 0-ball is a point and an oo -ball is an open half-space. An a-ball b is empty if b D S = 0. 
Any subset T C S of size \T\ = k + 1, with < k < 3, defines a k -simplex that is the convex 
hull of T, also denoted by conv(T). The general-position assumption assures that all A;-simplices 
are properly ^-dimensional. For < k < 2, a A;-simplex is said to be a-exposed if there is 
an empty a-ball b with T = db D S } where db is the sphere or plane bounding b. A fixed a thus 
defines sets Fk, a of a-exposed A;-simplices for < k < 2. The a-shape of S, denoted by <S„, is the 
polytope whose boundary consists of the triangles in i*2,co the edges in Fi jCi} and the vertices in 
Fo, a (see figures 1 and 2). The A;-simplices in Fk, a are also called the k-faces of S a . 

We still need to specify which connected components of 1R 3 — dS a are interior and which are 
exterior to S a . Fix the value of a and notice that for each a-exposed triangle there are two 
(not necessarily empty) a-balls, bi ^ b 2} so that T C dbi and T C db 2 . If both a-balls are empty 
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then (Tj does not belong to the boundary of the interior of S a . Otherwise, assume that bi is empty 
and that b 2 is not. In this case, bounds the interior of S a . More specifically, the interior of 
S a and the center of bi lie on different sides of aj. The definition of interior and exterior of S a is 
possibly more natural in the context of Delaunay triangulations and a-complexes as described in 
section 3. 
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Figure 1: Two tori, [n = 800, |F 2 | = 12197] 
The points are randomly generated on the surface of two linked tori. Six different a-shapes for 
values of a decreasing from top to bottom and left to right are shown. The first shape is the 
convex hull, for a = +oo; the last shape is the point set itself, for a = 0. The a-value used in the 
forth frame neatly separates the two tori. Further decreasing a disassembles the shape. Singular 
triangles, which do not bound the interior of the shape, are shown in darker color. 
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Figure 2: Bust, [n = 2630, |F 2 | = 35196] 
This point set is based on a demo data set for Silicon Graphics' Solidview program, of course, 
without any of the original connectivity information. The erasing sphere is shown to the right of 
the shape. Apart from a dense conglomerate of points representing part of the person's brain and 
brain stem, the set is basically hollow with most points representing skin. 



5 



3 Related Geometric Concepts 



There are quite a few natural geometric concepts that are closely related to a-shapes. Some of them 
are discussed in this section. In each case, the emphasis is on how the concept is related to a-shapes 
and how this relation can enrich our understanding of a-shapes. Section 3.1 discusses a- hulls and 
a-diagrams. Section 3.2 briefly reviews Delaunay triangulations and their dual incarnations known 
as Voronoi diagrams. The relevance of the Delaunay triangulation of a point set is that each a- 
shape of the set is the underlying space of a subcomplex of the triangulation. These subcomplexes 
are termed a-complexes and defined in section 3.3. Extensions of these notions are mentioned in 
section 3.4. 

3.1 Alpha Hulls and Alpha Diagrams 

Recall from section 2 that for positive real a an a-ball is defined as an open ball with radius a. 
For a = 0, it is a point, and for a = oo, it is an open half-space. Given a finite point set S C 1R 3 } 
an a-ball is empty if b D S = 0. For < a < oo, the a-hull of S, denoted by 7i a} is defined as 
the complement of the union of all empty a-balls. This is the formal counterpart of the styrofoam 
object described in section 2. Sample members of the continuous family of a-hulls are the convex 
hull of S, for a = oo, and S itself, for a sufficiently small. Observe that 7i ai C 7i a2 if a\ < ct 2 . 

Another interesting concept defined by a-balls is what we call the a-diagram of S, denoted by 
U a . For < a < oo, U a is the union of all a-balls whose centers are points in S 3 Observe that 
a point x £ 1R 3 belongs to U a iff the a-ball centered at x is not empty. Denote this a-ball by b x . 
This implies the following close relationship between 7i a and U a . 



Consider the boundary of U a . It consists of spherical caps, circular arcs, and vertices which we 
call corners. These are the 2-, 1-, and 0-faces of U a . These caps, arcs, and corners are in close 
correspondence with the vertices, edges, and triangles of S a . Some definitions are needed to 
describe this correspondence. 

Let a be fixed, with < a < oo, and let T be a subset of S of size \T\ = k + 1, with < k < 2. 
Define Kj = f] peT db p} where b p is the a-ball centered at p } as before. Besides general position 
assume Kj ^ 0. For \T\ = 1, Kj is a sphere; for \T\ = 2, Kj is a circle; and for \T\ = 3, Kj is a 
pair of points. It follows from the definitions that T is a-exposed iff Kj contains at least one face 
of the boundary of U a . If \T\ = 1 this face is a cap; if \T\ = 2 it is an arc; and if \T\ = 3 it is a 
corner. This fact can be expressed as follows. 

(Tj is a &-face of S a , with < k < 2 <^=> Kj contains at least one (2 — &)-face of U a . 

3 In chemistry and biology, a-diagrams are known as space-filling diagrams. However, there they are usually 
not restricted to equally large balls. This restriction can be removed with weighted a-shapes and a-diagrams (see 
section 3.4). 



x £ U a 
x £ "H, 



a 



■a 




and 



6 



Moreover, the number and structure of (2 — &)-faces contained in Kj are reflected by and the 
way it is embedded in S a . For example, if \T\ = 3, then Kj contains no, one, or two corners of 
U a . First, it contains no corner iff is not a triangle of S a . This is mentioned above. Second, 
Kt contains one corner iff bounds the interior of S a . Third, both points of Kj are corners of 
U a iff both sides of face the outside of <S„, that is, is a triangle of S a that does not bound 
its interior. Such a triangle will be called singular in section 5.2. Similarly, if \T\ = 2 and is an 
edge of <S„, then every angular interval between two incident triangles of S a that faces the outside 
corresponds to an arc of U a contained in Kj. If there is no incident triangle (in this case, is an 
singular edge) then the entire I-sphere Kj belongs to the boundary oiU a . Analogous statements 
can be made for \T\ = I where spherical caps on Kj correspond to solid angles around the vertex 
(Tj that face the outside of S a . 

All this amounts to a one-to-one correspondence between the (2 — &)-faces of U a and the &-faces 
of <S„, provided the latter are interpreted with multiplicities reflecting the number of exposed 
sides, angular intervals, or solid angles. The one-to-one correspondence also preserves incidences. 
This suggests that lA a , or more specifically, the boundary of U a , be represented by <S„, or more 
specifically, the faces of S a and their incidences. 

3.2 Delaunay Triangulations and Voronoi Diagrams 

A finite point set S C 1R 3 defines a special triangulation known as the Delaunay triangulation of S 
(see for example [11, 35]). Assuming general position of the points, this triangulation is unique and 
decomposes the convex hull of S into tetrahedra. The triangulation is named after the Russian 
geometer Boris Delaunay (also Delone) who introduced it in his seminal paper [6] in 1934. As 
explained below, the Delaunay triangulation of S is dual to another complex defined by S known 
as the Voronoi diagram [41,42]. Both complexes are related in an interesting way to the family 
of all a-shapes of S. The relationship between a-shapes and Delaunay triangulations will be of 
particular importance for this paper. 

Delaunay Triangulations. For < k < 3, let Fk be the set of A;-simplices aj = conv(T), T C S 
and \T\ = k + 1, for which there are empty open balls b with dbD S = T . Notice that F = S. The 
Delaunay triangulation of S } denoted by T> } is the simplicial complex defined by the tetrahedra 
in F 3} the triangles in F 2} the edges in and the vertices in F . By definition, for each simplex 
(Tj G £>, there exist values of a > so that (Tj is a-exposed. Conversely, every face of S a is a 
simplex of T>. This implies the following relationship between the Delaunay triangulation and the 
boundary of S a . 

F k = U F K«, for < & < 2. 

0<a<oo 

We take advantage of this relationship by representing the family of a-shapes of S implicitly by 
the Delaunay triangulation of S. This will be described in detail in section 5. 

Voronoi Diagrams. For a point p G S } define V(p), the Voronoi cell of p } as the set of points 
x G M 3 so that the Euclidean distance between x and p is less than or equal to the distance 
between x and any other point of S. Each Voronoi cell is a convex polyhedron, and the collection 
of all Voronoi cells, one for each point of S } defines the Voronoi diagram of S } denoted by V. We 
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call a Voronoi cell also a 3 -cell of V. Each 2 -cell of V is the intersection of two Voronoi cells; 
each 1-cell or edge is the intersection of three 3-cells; and each 0-cell or vertex is the intersection 
of four 3-cells. There is a natural one-to-one correspondence between the ^-simplices of T> and 
the (3 — &)-cells of V. Let T be a subset of S of size \T\ = k + 1, with < k < 3, and define 
V T = f] P e T V(p). 

(Tj is a £;-simplex of T> <^=> Vj is a (3 — A;) -cell of V, for < k < 3. 

This correspondence preserves (or reverses) incidences which implies that T> and V are indeed dual 
to each other. 

Observe that Vj is the set of all points x for which there exists an empty open ball b x centered 
at x with T C db x D S. Equality holds iff x belongs to the relative interior of Vj. It follows that 
(Tj is a-exposed iff there is a point x in the relative interior of Vj whose distance from the points 
in T is a. Since Vj is convex there is a single interval so that is a-exposed iff a belongs to 
this interval. We will exploit this fact later when we discuss how to decide when a simplex of T> 
is a-exposed. 

3.3 Alpha Complexes 

Since all faces of S a are simplices of T> } it follows that the interior of S a is naturally triangulated 
by the tetrahedra of T>. This idea leads to the concept of a-complexes as defined shortly. A 
(three-dimensional) simplicial complex is a collection C of closed ^-simplices, for < k < 3, that 
satisfies the following two properties. 

(i) If (Tj G C then (Tji G C for every T' C T . In other words, with every simplex <tx, C contains 
all faces of (Tj as well. 

(ii) If <tx, (Tji G C, then either (Tj D (Tj 1 = 0, or (Tj D (Tj 1 = CTnT' = conv(T fl T'). Note that (i) 
implies that this face is also in C. In other words, the intersection of any two simplices in C 
is either empty or a face of both. 

A subset C C C is a subcomplex of C if it is also a simplicial complex. 

Each £;-simplex (Tj of T> defines an open ball bj bounded by the smallest sphere dbj that contains 
all points of T . Let qt be the radius of bj. For k = 3, dbj is the circumsphere of cr^; for k = 2, the 
circumcircle of (Tj is a great circle of dbj] and for k = 1, the two points in T are antipodal on dbj. 
Call dbx the smallest circumsphere and qt the radius of <tx- For 1 < k < 3 and < a < oo, define 
Crjfe )Q: as the set of A;-simplices (Tj G © for which &x is empty and < ct. Furthermore, define 
Co, a = S } for all a. The sets Gk, a do not necessarily define a simplicial complex because it can 
happen that Gz, a contains a tetrahedron but not all triangles of this tetrahedron belong to G^,a- 
Similarly for triangles and edges. With this in mind, we define the a-complex of S } denoted by 
C a , as the simplicial complex whose ^-simplices are either in Gk, a or they bound (k + l)-simplices 
of C a . By definition, C ai is a subcomplex of C a2 if ct\ < ct 2 . 

The underlying space of C a , denoted by \C a \, is the union of all simplices of C a , or in other words, 
the part of 1R 3 covered by C a . Thus, the underlying space of C a is a polytope in the sense specified 
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in section 2. Indeed, we have the following most important property of C a , which we present 
without proof. 

For all < a < oo, S a = \C a \. 

This can be considered an alternative definition of a-shapes. It makes it easy to specify the intervals 
of the simplices of T> alluded to in the above discussion of Voronoi diagrams. For example, let 
be a £;-simplex of T>. If bj is empty then belongs to C a iff a £ (qt, oo]. 4 

3.4 Extensions 

The definitions presented in sections 2 and 3 above can be extended in various ways. So far we 
refrained from mentioning these extensions in order to avoid unnecessary complications and to be 
faithful to the currently available implementation of the concepts in this paper. For completeness, 
negative values of a, weighted points, and higher dimensions are now briefly discussed. 

Negative Alpha Values. This extension has been described in [14] for the two-dimensional 
case. 5 For negative a, a-complexes are most naturally defined as subcomplexes of the so-called 
furthest-point Delaunay triangulation of S (see for example [11,35]). For T C S and \T\ = 4, 
the tetrahedron belongs to this triangulation iff bj contains all points of S — T . The a- 
shape is, again, the underlying space of the a-complex. These shapes exhibit far less interesting 
geometric and topological properties than the ones for positive a, and are thus less interesting for 
applications. We omit further details. 

Weighted Points. Recall the relationship between a-shapes and a-diagrams described in sec- 
tion 3.1. It is interesting to consider diagrams for different ball sizes, and this is indeed done in 
chemistry and biology where space-filling diagrams are usually defined as unions of balls with ar- 
bitrary and thus possibly different radii. In order to represent such diagrams by polytopes similar 
to a-shapes it is necessary to introduce weighted a-shapes. These can be defined using subcom- 
plexes of so-called regular triangulations (see for example [12,28]). Given a finite set of points, 
each with a real weight, the regular triangulation is a unique simplicial complex whose underlying 
space is the convex hull of the point set. If all weights are the same then it equals the Delaunay 
triangulation of the points. Details can be found in [13]. 

Higher Dimensions. It is fairly straightforward to generalize all important concepts of this 
section (like a-shapes, a-hulls, a-diagrams, a-complexes, Delaunay triangulations, and Voronoi 
diagrams) to finite point sets S in M d , for arbitrary dimension d. This generalization, combined 
with an extension to weighted points is developed in [13]. Note however, that the implementation 
details become progressively "hairier" with increasing dimension, and the worst-case complexity 
of the problem grows exponentially. For example, if S is a set of n points in M d then the Delaunay 
triangulation of S can consist of up to ©(n^ - ^-!) faces (see [38]). Although the running time of 
the programs constructing a-shapes will get substantially worse as d increases, there might be 
applications that warrant implementations in low dimensions higher than three. 

4 Because of the general-position assumption we have a ^ qt ■ It is therefore irrelevant whether the interval is 
open or closed at its left endpoint. 

5 The definitions in [14] are slightly different from the ones in this paper. In particular, their a is the same as 
minus one over a in this paper. Thus, our negative values of a correspond to their positive values. 
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4 Combinatorial Analysis 



In contrast to a-hulls and a-diagrams, the a-shapes of a finite point set form a discrete family, 
even though they are defined for all real numbers a, with < a < oo. Indeed, S ai ^ S a2 iff 
Ufc=o Gk, ai 7^ Ufc=o Cfc,a 2 - Thus, S ai 7^ S a2 iff there is an empty open ball bounded by a smallest 
circumsphere of an edge, triangle, or tetrahedron of T> whose radius lies between a 4 and a 2 . Such 
a radius is referred to as an a-threshold because it separates two a-shapes. The number of a- 
shapes exceeds the number of a-thresholds by one. It follows that one plus the total number of 
^-simplices of T> } for 1 < k < 3, is an upper bound on the number of different a-shapes. An upper 
bound on the number of simplices of T> can be obtained using a relationship between Delaunay 
triangulations in 1R 3 and certain convex polytopes in 1R 4 . We briefly describe this relationship in 
the next paragraph and refer to [11,38] for details. 

Lifting Map. Identify 1R 3 with the Xix 2 x 3 -sp&ce in 1R 4 } that is, the subspace x 4 = 0. The lifting 
map is a geometric transform that projects points p = (7Ti, 7r 2 , 7r 3 ) in 1R 3 along the x 4 -axis onto 
the paraboloid of revolution U: x 4 = Y^=i x 1 m ^? 4 - Let pu = (tti, 7r 2 , 7r 3 , Y1^=i ^J) be the image 
of p } and define Su = {pu \ P £ S}. The convex hull of Su, conv(Su), is a convex polytope with n 
vertices in M 4 . A facet belongs to the lower boundary of this polytope if the polytope lies on the 
side of the positive x 4 -axis of the hyperplane that contains the facet. Otherwise, it belongs to the 
upper boundary of conv(Su)- If we project all facets of the lower boundary of conv(Su) parallel 
to the x 4 -axis into 1R 3 } along with their subfaces, then we obtain the Delaunay triangulation of S 
(see [11]). 

Upper Bounds. According to the upper bound theorem for convex polytopes, the maximum 
numbers of 1-, 2-, and 3-faces of a convex polytope with n > 5 vertices in M 4 are \{n 2 — n), n 2 — 3n, 
and \{n 2 — 3n), respectively (see for example [3,32]). The lifting map implies the same upper 
bounds for with 1 < k < 3. As a matter of fact, the upper bound for \F 3 \ is one less than 

for the number of 3-faces of conv(Su) because at least one 3-face belongs to the upper boundary 
of conv(Su)- By a result of [38], these bounds are tight even though the vertices of conv(Su) are 
constrained to lie on a second-degree surface in M 4 . We summarize these results for n = \S\. 

\Fo\ = n, |Fi| < i(n 2 - n), \F 2 \ < n 2 - 3n, and |F 3 | < ^(n 2 - 3n - 2). 

By adding one to the sum of the bounds for \F 2 \ } and I-F3I, we obtain the following result. 

S has at most 2n 2 — 5n different a-shapes. 

This bound is too pessimistic for two reasons. First, although the upper bounds on the number of 
simplices of T> are tight, there are many fewer simplices for most point sets. For example, if the 
n points are uniformly distributed in the unit ball then the expected number of simplices is only 
0(n) (see Dwyer [9]). Second, not all edges and triangles of T> have a smallest circumsphere that 
bounds an empty open ball. However, since the circumsphere of every tetrahedron of T> bounds 
an empty open ball by definition, the number of different a-shapes is always at least a fraction of 
the number of simplices of T>. 
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5 Algorithms 



As described in section 3, the family of a-shapes of a finite point set S can be represented by 
the Delaunay triangulation of S. In this representation, each simplex of T> is associated with 
an interval that specifies for which values of a the simplex belongs to the a-shape. Section 5.1 
discusses the construction of T> } and section 5.2 explains how the intervals of the simplices are 
computed. For completeness, section 5.3 gives the formulas that can be used to implement the 
required primitive operations. 

5.1 Three-dimensional Delaunay Triangulations 

The construction of Delaunay triangulations is a popular topic in the area of geometric algorithms 
(see for example [11,35]). Indeed, various different approaches have been studied and described 
in the literature. Some approaches are based on the lifting map mentioned in section 4, which 
transforms the problem into one of constructing the convex hull of a four-dimensional point set. 

The algorithm adopted for our implementation of a-shapes has been suggested by Joe [25]. It is 
based on the idea of local transformations or flips. 6 The algorithm can be viewed as a generalization 
of the edge-flip method for two-dimensional triangulations by Lawson [27]. Note, however, that 
the straightforward generalization of the two-dimensional algorithm to 1R 3 fails to always compute 
the Delaunay triangulation (see [24]). Nevertheless, the correctness of the flip algorithm in 1R 3 
can be established if the points are added one by one. We sketch the structure of the resulting 
incremental algorithm below and discuss the notion of a flip later. 

Incremental-flip Algorithm. 

Sort S along some fixed direction, 

and relabel the points so that (pi 7 p 2} . . . ,p n ) is the sorted sequence. 

Initialize T> to the triangulation whose only tetrahedron is <tx, T = {piiP2iPziP±} ■ 

for i := 5 to n do 

Add pi by connecting it to all vertices, edges, and triangles of T> visible from pi. 
Use flips to transform T> to the Delaunay triangulation of {pi 7 p 2} . . . 7 Pi}- 
end for. 

Consider two tetrahedra (Tji and (Tjh that share a common triangle aj. For example, T = 
{pi } p 1} pk} } T' = T U {p u }i and T" = T U {p v }. Triangle is called locally Delaunay if it 
belongs to the Delaunay triangulation of T' U T" . This is the case iff the point p v lies outside 
the sphere dbji. Local Delaunayhood is a necessary condition for (Tj to belong to the Delaunay 
triangulation, but it is not sufficient. Nevertheless, Delaunay [6] proved that if all triangles are 
locally Delaunay then the triangulation is a Delaunay triangulation. 

The incremental-flip algorithm needs to restore Delaunayhood whenever a new point is added 
to the triangulation. For this, it identifies triangles that are not locally Delaunay and tries to 
flip them. Let (Tj be such a triangle bounding the tetrahedra (Tji and (Tjn. Again, assume 
T = {pi } p 1} pk} } T' = T U {p u }i and T" = TU {p v }- If cx' U (Tt" is convex then the triangle (Tj can 

6 As recently shown, the flip algorithm can be extended to compute regular triangulations in arbitrary dimensions 
[16]. Regular triangulations are useful for weighted a-shapes (see section 3.4 or [13]). 
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be replaced by the edge connecting p u and p v . Together with this edge, the triangles connecting it 
with the three vertices of T are added. This operation is called a triangle-to-edge flip. Otherwise, 
(Tji U (Tt" is not convex. Assume there is a third tetrahedron <7x»' that is spanned by four of the 
five points of T' U T", for example, T'" = {pi } p 17 p U} p v }. In this case, the three triangles incident 
to the edge connecting pi and pj can be replaced by a single triangle with endpoints p U} p V} and 
Pk- We call this an edge-to-triangle flip. If <7x»' does not exist then there is no flip that can remove 

<7y. 

In spite of possible triangles that are neither locally Delaunay nor can be flipped according to the 
above rules, the correctness of the increment al-flip algorithm can still be established [25]. In the 
worst case, it takes time and storage 0(n 2 ) } where n = \S\. Experiments provide evidence that 
it performs significantly better for most point sets. However, the worst case of 0(n 2 ) cannot be 
avoided because T> can have up to a quadratic number of simplices (see section 4). 



5.2 Intervals and Face Classification 

For each simplex Of £ T> there is a single interval so that Of is a face of the a-shape S a iff a is 
contained in this interval. This was mentioned in section 3. It will be convenient to study these 
intervals for the a-complex C a rather than the a-shape. Also, we break each interval into three 
(possibly empty) parts that correspond to values of a for which the simplex is an interior, regular, 
or singular simplex of C a . 

A simplex Of £ C a is said to be 
interior if Of ^ <9<S„, 

regular if Of £ dS a and it bounds some higher-dimensional simplex in C a , and 
singular if Of £ dS a and it does not bound any higher-dimensional simplex in C a . 

Notice that there are Delaunay edges and triangles that can never be singular because their 
smallest circumsphere encloses other points of S. 7 Therefore, we call a simplex aj £ T> 

{attached if \T\ = 2, 3 and bj H S ^ 0, and 
unattached otherwise. 

Recall that qt is the radius of the smallest circumsphere of Of. In order to break up the interval 
for which Of belongs to C a , we introduce values // and JI T for which Of changes from singular to 
regular and from regular to interior, respectively. Before that, however, let up(o"x) be the set of 
all simplices in T> that contain a simplex Of £ T> } with \T\ < 3, as a proper face, that is, 

up(<r T ) = {(T T , £ V | T C T'}. 

7 It is convenient to extend the general-position assumption so that no smallest circumsphere of two or three 
points of S contains another point of S. A slightly more general assumption is the following. If a sphere contains 
three points of S then no two of them are antipodal, and if it contains four points then no three lie on a great-circle 
of the sphere. 
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singular 


regular 


interior 


tetrahedron 






( ot. oo] 


edge or triangle, 4_ dconv(S'), unattached 
4 9conv(S'), attached 
G dconv(S'), unattached 
G dconv(S'), attached 


(qt,h t ) 

(Qt,H t ) 


(//_,, JJ't ) 

(n T ,°°\ 


(/I T ,oo] 
(ji T , oo] 


vertex, 4_ dconv(S) 
G dconv(S') 


[0,/£ T ) 




(/I T ,oo] 



Table 1: Intervals of a values for which G D belongs to the a-complex C t 



If (Tj is a tetrahedron, define fi T = fi T = qt- Otherwise, 

fi = minj^x' | ct' £ upl^x), unattached} and 
JI T = maxj^x' | cx' G up(o"x)}. 

It is sufficient to consider only the set 

u Pi(°"x) = {o"x' £ up(cr T ) | |T'| = |T| + 1}, 

that is, all simplices incident to whose dimension is one higher than that of <tx, in order to 
derive the values // and ~p T ; 

■ ( {qt 1 I cx' £ uPil^x), unattached} \ 
— T mm y U {fJ- T , | cx' £ uPil^x), attached} J ' 

and 

/Z T = max{/Z T ; | G up 1 (o"x)}. 

Specifying Intervals. The intervals of a values in which is an interior, regular, or singular 
simplex of C a are shown in the table 1. Because of the general-position assumption, a is different 
from all g values and therefore also from all // and JI values. We can thus define all intervals as 
open, except at endpoints and oo. It is necessary to distinguish simplices that bound the convex 
hull of S from the others. The next paragraph briefly explains the entries of table 1 for the case 
of triangles that do not bound the convex hull of S. The arguments for tetrahedra, triangles on 
the convex hull, edges, and vertices are similar. 

Consider a triangle G T> } T = {pi } p 1} pk} } that does not bound the convex hull of S; we 
denote this by 4_ dconv(S'). Let (Tji and (Tjh be the two incident tetrahedra in T> } and assume 
T' = T U {p u }i and T" = T U {p v }. Furthermore, let < qt> < Qt" < oo; in other words, // = qji 
and ~p T = g T n. Now, fix a value for a. If qt" < a < oo then the triangle is not a-exposed. 
It will, however, be part of the interior of <S„, because both incident tetrahedra are in C a . If 
qt' < a < qt" then the triangle is a-exposed and (Tji is in C a but (Tjh is not. This means that (Tji 
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is a regular triangle of C a . For a < qt* } neither (Tji nor (Tjh are tetrahedra of <S„, but can still 
be a singular triangle, that is, iff qt < a and neither p u nor p v are inside bj. If one of the two 
points is inside &x, then is attached, and can never be a singular triangle of C a , no matter 
what a value is selected. 

The a-complex consists of all interior, regular, and singular simplices for a given a value. The 
interior of the a-shape is triangulated by the interior simplices. The boundary of the interior is 
formed by the set of regular triangles and their edges and vertices. 

Consistent with the definition in section 4, we refer to the endpoints of the intervals in table 1 
as a-thresholds. This does not include and oo. Since all // and JI T values are g values of 
other simplices, each a-threshold is the radius of a simplex in T>. More specifically, the set of 
a-thresholds is exactly the set of radii of all unattached ^-simplices for 1 < k < 3. Define the 
a-spectrum as the sorted sequence of a-thresholds. This concept will appear again in section 6. 

Computing Intervals. Assume that each simplex £ T> is marked as either "£ dconv(S')" 
or dconv(S')" after the construction of T>. With this, the above intervals can be computed by 
classifying as attached or unattached, and by computing qt, // , and /Z T , whenever applicable. 
We said that is attached iff one of the simplices that contain has a vertex in &x, the 
open ball bounded by the smallest circumsphere of aj. This implies that can be classified in 
time proportional to |up 1 (o"x)|. The time it takes to classify all simplices is proportional to the 
number of simplices in T> } because each simplex has only a constant number of faces. In other 
words, assuming that constant time suffices to decide whether or not a point belongs to bj (see 
section 5.3), a simplex can be classified in constant amortized time. 

Furthermore, assume that, given T with £ T> } qt can be computed in constant time (again, see 
section 5.3). By processing tetrahedra before triangles before edges before vertices, we can get // 
and ~p T simply as the minimum and maximum of the values qt' } /« t ,, and /Z T; , for (Tji £ up 1 (o"x). 
This also takes only constant amortized time per simplex. 

5.3 Geometric Primitives 

What are the primitive operations needed to compute a-shapes in _ff? 3 ? Constructing Delaunay 
triangulations requires two geometric tests. These are a test for deciding on which side of a plane 
spanned by three points a fourth point lies, and one for deciding on which side of a sphere spanned 
by four points a fifth point lies. In order to generate the intervals of table 1, we need to compute 
the radius of the smallest circumsphere of a tetrahedron, triangle, or edge, and test whether a point 
lies inside or outside this sphere. While the two tests required for Delaunay triangulations are 
fairly common in geometric algorithms and computer graphics, the operations involving smallest 
circumspheres of triangles and edges are rather specialized. All operations share the problem 
of degenerate cases, which we can ignore because of the general-position assumption (see also 
section 6). This section gives a formula for each of the primitive operations mentioned above. 

Assume that the points of S are labeled as pi 7 p 2} ■ ■ ■ } p n} and that each point pi is given by the 
vector (vr 8 ,i, 7T 8 ' j2 , vr 8 , 3 ) of its three coordinates. To simplify the notation in the remainder of this 
section, we use minors, which are determinants of submatrices of a given matrix. For convenience, 
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define 7r 8j0 = 1 for all i and use the following notation for minors: 



A/f 1 ' 8 . 2 '-' 8 * = det 



/ "'liji 7r «lj2 • • • "''bit \ 



V X 'tJl 7F UJ2 • • • ^ikJk / 



Indeed, all geometric primitives are expressed in terms of determinants, which provide a convenient 
and compact notation. This does not exclude the possibility of implementing the primitives using 
equivalent formulas; these would be, in some sense, ways of evaluating determinants that differ 
from the usual constructive definition. 

Plane Test. Let T = {pi } p 17 Pk} and define hj as the unique plane that contains all three 
points of T . This plane can be oriented if we replace the set T by the sequence T; for example, 
T = (pi } p 1} pk). Then one side (or open half-space) of hj can be called positive and the other 
negative. We also refer to these as the positive and negative sides of the sequence (pi } p 1} pk). 

T = (pi, pj , pk): p u lies on the positive side of hj <^=> •M-i^so > 0? (5-1) 

and p u lies on the negative side if the determinant is negative. Intuitively, p u sees the sequence of 
three points pi } p 1} pk in a clockwise order iff p u lies on their positive side. Similarly, p u sees the 
sequence in a counterclockwise order iff p u lies on the negative side of the points. The sign of the 
determinant is called the orientation of the sequence (pi } p 1} pk } p u )- Notice that the determinant 
equals zero iff the points are in degenerate position, that is, they lie on a common plane. Observe 
also that the orientation of a permutation of a sequence of four points is the same as the orientation 
of the sequence itself, provided the number of transpositions is even. Otherwise, it is the opposite. 
This follows trivially from the fact that the value of the determinant changes sign whenever two 
rows are exchanged. 

Sphere Test. Given a set T = {pi } p 1} pk } p u } } we need to decide whether another point p v lies 
inside or outside the sphere dbj. We can assume that the degenerate case where p v £ dbx does 
not occur. A possible implementation of this test is discussed in [15] using an extension of the 
lifting map (as mentioned in section 4) to three-dimensional spheres. Consider the paraboloid 
of revolution U: x 4 = J2i=i x i i n ^? 4 - A sphere db with center c = (71,72,73) and radius p is 
mapped to the hyperplane dbu'- x 4 = Y3_=i(^li x i ~ 7?) + P 2 ■ This hyperplane has the property 
that U H dbu projected along the x 4 -axis into the Xix 2 x 3 -sp&ce yields db. Moreover, a point p lies 
inside (outside) db iff pu lies vertically below (above) dbjj. The resulting formula assumes that 
each point pi £ S has a fourth coordinate 7r 8 ' j4 = Yfj=i ^1 3 - 

T = {Pi,Pj,Pk,Pu}- Pv lies inside db T <i==> M\fy}) ■ M\fy^ > (5-2) 

The first minor, M.^^ioi 1S a corrective term that is necessary because the sphere does not change 
if the first four points are permuted. The second minor, ^(^34 v , expresses the fact that the 
lifting map transforms a sphere test in 1R 3 to a hyperplane test in JR A . 

Radius of a Smallest Circumsphere. Next, we consider computing the radius qt of dbx-, the 
smallest circumsphere of <tx, for all A;-simplices Of £ T> } with 1 < k < 3. The formulas for the 
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square of qt are given in (5-3) through (5-5). Note that computing q\ will be sufficient for our 
purposes since qt can never be negative. We distinguish the cases when k = \T\ — 1 is 1, 2, or 3. 



T = {pi, Pj }: e 2 T = ~ ~ " — (5-3) 



T = {pi, Pj , Pk }: g 2 T = 7 v . y 2 . v (5-4) 

4-((^ 2 ',to) +(m\$ ) +(m\% ) ) 

.i,i,fc,-u\ 2 . / ,.i,j,k,u\ 2 . ( ,.i,j,k,u\ 2 . , K A,j,k,u i.i,j,k,u 



/ M i,3,k,u \ I M i, 3 ,k,u \ I M i, 3 ,k,u \ . KAi,],k,u p.%,],k,u 

m f i2 I 2,3,4,0i ^IV^' 1,3,4,0,1 T I 1,2,4,0 i + 41 -^'l, 2,3,0 -^'l, 2,3,4 

T = {p t ,Pj,Pk,Pu}- Q T = ; — — —2 (5-5) 

4 • (M\%$ 



In order to explain these three formulas we introduce some notation. Let a, 6, c, J be points in 
M 3 . We write |a&| for the length of the edge conv({a,6}), and \abc\ for the area of the triangle 
conv({a, 6, c}). In the case of two points, qt is the same as half the distance between pi and p r 
Equation (5-3) follows because 

\piPi\ 2 = {M^ y+ (M^ y+ (M^ y. 

To handle the case of a triangle, that is, T = {pi } p 1} pk} } we use the formulas 



\PiPj I ' \PiPk I • \PkPi I , 
g T = j j and 

4 • \PiP 3 Pk | 

\piPiPk\ 2 = \ ■ ((^2',? ) 2 + {M^ y + (m^ 

which can be found in any good mathematical handbook. Finally, we obtain (5-5) using the 
extension of the lifting map mentioned above. If db is the sphere through points pi } p 1} pk } p u then 
dbjj is the hyperplane through the four points P i,u , Pj,u , Pk,u , Pu,u • The equation for the hyperplane 
can be computed directly from the coordinates of the lifted points. From this equation it is easy 
to compute the center and the radius of db. 

Attached and Unattached Edges and Triangles. We still have to consider the problem of 
deciding whether an edge or triangle (Tj £ T> is attached or not. By definition, is attached if 
there is a <tr £ upi(crx) so that the point in R — T belongs to bj. If is an edge, say, T = { P i,Pj} 
and R — T = {pk}, this can be done by comparing qt with the distance between pk and Pl+ ^ . 
Straightforward algebraic manipulations lead to the following equation. 

T={p i ,p j }: Pk eb T ^ f:{M^y-f:(M^+M^y>o (5-6) 

£=1 £=1 

Now let (Tj be a triangle, for example, T = {pi 7 p 17 Pk} and R — T = {p u }- To see whether or not 
the point p u belongs to &x, we compute the center c of the circumsphere dbn of the tetrahedron 
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<tr. Observe that p u G ^ iff c and p u do not lie on the same side ol the plane through pi } p 1} pk. 
In other words, we need to test whether or not the sequences (p U} pi } p 1} pk) and (c } pi } p 1} pk) 
have different orientation. Some rather tedious algebraic manipulations are needed to derive the 
following equation which expresses the derivation in terms of minors. 

{KAi,3,k,u XAhhk , KAi,3,k,u XAhhk 
yW 2,3,4,0 JVl 2,3,0 T • /w 1,3,4,0 JVl l,3,0 
(5-7) 
M i,j,k,u M i,j,k _ r, _ M i,j,k,u M i,j,k r, 
' JV1 1, 2,4,0 JV1 1,2,0 L JV1 1, 2,3,0 yvl l,2,3 ^ u 

General Position Revisited. The general-position assumption used in this paper assures that 
no geometric test is ambiguous. We summarize and revise the necessary assumptions below and 
include pointers to the formulas for which the assumptions are relevant. 

• No 4 points lie on a common plane; compare with (5-1). 

• No 5 points lie on a common sphere; compare with (5-2). 

• No smallest circumsphere of 2, 3, or 4 points has a radius equal to any given a; compare 
with (5-3), (5-4), and (5-5). 

• No point lies on the smallest circumsphere of 2 or 3 other points; compare with (5-6) 
and (5-7). 

These assumptions are indeed very restrictive and rarely true for real-life data. We will deal with 
this apparent shortcoming in section 6.2. 
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6 Implementation 



Our current implementation ol a software tool for a-shapes in 1R 3 consists of the following three 
parts. 

1. A program that constructs Delaunay triangulations using flips (see section 5.1). 

2. A program that computes the a-intervals for all simplices in a Delaunay triangulation, and 
then sorts the endpoints of these intervals (see section 5.2). 

3. An a-shape visualizer that enables the user to manually select different a values and render 
the corresponding shape on a graphics workstation (see figures 1 through 5). 

Parts 1 and 2 are preprocessing steps that take time 0(n 2 ) and 0(ralog m), where n is the number 
of points and m is the number of simplices of T>. The current code for part 3 takes time 0(m) 
to render a particular a-shape. Improvements based on fast data structures for intervals are 
forthcoming. 

One important aspect of the implementation is its robustness. By this we mean that the program 
will either produce the correct output for the original data set S, or, in case of degeneracies, it 
guarantees to give the correct result for a set S(e) arbitrarily close to the original input. This is 
achieved by a symbolic perturbation scheme briefly described in section 6.2. The method avoids 
possible conflicts between the topological and geometric structure of the data by using exact 
arithmetic and by perturbing the original data set such that all degeneracies disappear. The 
program can thus be considered to be purely combinatorial and logical so that correctness in the 
strict sense is possible in principle. Note that the perturbation is infinitesimal and symbolic, that 
is, S and S(e) can be viewed as arbitrarily close together, and the computational overhead is 
independent of e. 

6.1 Data Structures 

There are two main data structures needed for storing the family of a-shapes of a given data set. 
One represents the connectivity and order among the simplices of the three-dimensional Delaunay 
triangulations. The other is used for the intervals assigned to the simplices of T>. A triangle-based 
data structure is used to store T>. This is briefly described in the paragraphs below. An interval 
tree can be used to store the collection of intervals (see for example [35]). The current version of 
our program, however, stores the a-spectrum using only a linear array. Recall that the a-spectrum 
is the sorted sequence of a-thresholds, and for practical reasons, the implementation adds and 
oo to this list. Universal hashing (see for example [5]) provides the link between the triangle 
structure and the array. 

The data structure used to store the three-dimensional triangulation is a specialized version of 
the edge-facet structure introduced by Dobkin and Laszlo [7]. Related data structures are the 
quad-edge structure due to Guibas and Stolfi [22], which can be used to model two-dimensional 
manifolds, and the cell-tuple structure by Brisson [2], which works in arbitrary dimensions. The 
edge-facet structure is designed for general cell complexes in three dimensions. By reducing the 
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scope to simplicial complexes, it is possible to improve the compactness and the speed of the 
structure. We refer to the result as the triangle- edge structure. 

The atomic unit of the triangle-edge structure is the so-called triangle- edge pair a = (cr,z), with 
< i < 5. It identifies six versions of the triangle <r, one for each of its six directed edges. 
Each triangle defines two edge rings. One edge ring traverses the edges of a in a counterclockwise 
order, the other traverses them in a clockwise order. Similarly, each edge defines two triangle 
rings traversing the incident triangles in the two opposite orders. Each triangle-edge a belongs to 
exactly one edge ring and exactly one triangle ring. 

The internal representation of the structure takes advantage of the fact that all edge rings have 
length three. It is thus possible to avoid actual pointers for the edge rings by merging the six 
triangle-edge pairs of two opposite edge rings into one record. Such a record allocates 30 (36) 
bytes per triangle, assuming that 2-byte (4-byte) integers are used as indices to the vertices, and 
4-byte integers for triangle-edge pairs. Further details are omitted. 

6.2 Simulated Perturbation 

For implementation purposes it is no longer appropriate to assume that the input points are in 
general position. This assumption would be too restrictive. On the other hand, in the context of 
three-dimensional a-shapes, it would be rather tedious to deal with the large number of special 
cases in ad hoc manner. For this reason, we apply a general technique, known as Simulation of Sim- 
plicity, or SoS [15], which acts as a black box between the implementation of a geometric algorithm 
and the input data. It allows a systematic treatment of all special cases on the level of geometric 
primitive operations. The SoS library consists of a set of carefully implemented primitives. It 
provides the programmer with the illusion of simple data while the actual input is in arbitrary 
and thus possibly degenerate position. 8 This section can only sketch the basic idea of SoS (refer 
to [15] for details). 

The idea of SoS is to perturb the given objects ever so slightly, in a manner that all degeneracies 
disappear. The perturbation should be small enough so that the nondegenerate position of objects 
relative to each other remains unchanged. Since it is usually rather difficult and costly to actually 
come up with such a perturbation, SoS performs it symbolically by substituting polynomials in e 
for the parameters specifying each object. In the context of this paper, the coordinate 7r 8J of the 
input point pi £ S is replaced by its perturbed version 

7r 8J (e) = nij + e(i,j), 

where e(i,j) = s 6 ** 3 and 8 is a sufficiently large constant. The choice of e(i,j) implies that general 
position of the perturbed input set S(e) is guaranteed if e is positive and arbitrarily small. Thus, 
e can be treated as an indeterminant in the symbolic evaluation of the primitive operations which 
are now based on S(e). 

The expressions that correspond to primitive operations are polynomials in e. The coefficients of 
these polynomials are expressions in the coordinates of the original data. For example, in order to 

8 The terms "simple," "general position," and "nondegenerate position" are used as synonyms. Notice how the 
"general case" of the algorithm designer is usually the "simple case" for the implementing programmer. 
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determine whether the sequence (pi } p 1} pk } p u ) } with i < j < k < u } has positive orientation, the 
sign of the polynomial 

M[f£ + M^l • 3) - ■ e(i, 2) + ■ e(i, 1) + ■ ■ ■ 

• • • + Mi& ■ e(k, 3)e(i, 1) + M$ ■ e(k, 3)e(j, 2) + e(k, 3)e(j, 2)e(i, 1) + • • • 

must be evaluated. Assume that the terms in the polynomial are sorted in order of increasing 
exponents of e. We say the evaluation has depth t if the coefficient of the it + l)-st term is the first 
one that does not vanish. Because e is assumed to be arbitrarily small, this term decides the sign 
of the entire polynomial. Notice, that the coefficient at depth is the same as the expression of the 
primitive for the original data. This implies that S(e) is consistent with S as far as nondegenerate 
configurations are concerned. Observe also that the polynomial never evaluates to zero since there 
is always a term with non-zero coefficient. In other words, S(e) is simple. 

The code that implements the polynomial of a given primitive operation can be generated auto- 
matically. The overhead for the symbolic perturbation as such is thus neglectable. However, the 
perturbation can only be simulated when exact arithmetic is used to compute the coefficients of 
the polynomials. More precisely, we need to be able to tell when a coefficient vanishes. Of course, 
exact arithmetic entails a somewhat higher cost for the low-level computations, but we believe that 
this is an adequate price for a compact and robust implementation of a rather involved geometric 
algorithm. 

Note that SoS does not allow the user to selectively remove certain types of degeneracies, while 
others remain. Rather, it fully implements what theoreticians call "general position." Observe 
that algorithms employing SoS produce results for the perturbed data, and not for the original 
one. Some postprocessing can be used to remove or repair parts of the output that collapse to a 
degenerate state because of the degeneracy of the input data. This is mentioned in [15]. As an 
example, consider the case of Delaunay triangulations and coplanar points on the boundary of the 
convex hull. The perturbation will move some points towards inside and some towards outside. In 
the Delaunay triangulation, the points that moved from the convex hull boundary inside will be 
covered by flat tetrahedra. These tetrahedra have infinitesimal thickness in the perturbed setting, 
and are completely flat in the original setting. It is easy to identify and remove these tetrahedra 
in a postprocessing step. 

6.3 Performance 

Table 2 shows the performance of the incremental-flip algorithm (see section 5.1) for a number 
of test runs. We count the number of flips, and the number of necessary evaluations of 5-by- 
5 and 4-by-4 determinants. The "max" and "mean" depth columns count the "depths" of the 
corresponding evaluations (see section 6.2). These columns give a measure of how degenerated 
a data set is. The experiments were run on Silicon Graphics workstations with 50 MHz MIPS 
R4000 CPUs and 48 Mb or more of main memory. Sample frames for the data sets molecule, 
tori, universe, bust, and phone_l can be found in figures 1 through 5. 

The running time for the presented examples is way better than predicted by the quadratic worst 
case. This is due to the fact that Delaunay triangulations usually do not reach their worst-case 
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upper bounds of 0(n 2 ) faces (see section 4). Table 2 shows running times that seem to be roughly 
proportional to n(logn) 2 . An exception to this rule are volumetric data with points on a regular 
grid or parallel slices (as in the data sets gridl, ku2, rat_T2a, and hsr). The incremental-flip 
algorithm is certainly less than ideal for such distributions and should possibly be replaced by a 
randomized version (see [16]). 9 

Observe the positive correlation between the number of determinant evaluations or long-integer 
operations and CPU seconds. Indeed, run-time profiles of the C code suggest that the majority 
of CPU cycles is used for the long-integer arithmetic computing determinants. We estimate that 
approximately 75% of the time is spent on long-integer arithmetic. The multiplication routine is 
responsible for more than half of it. 

Table 3 shows the performance of the program that generates the a-intervals for all simplices in 
the Delaunay triangulation (see section 5.2). The number of simplices in a triangulation is denoted 
by \T>\ = J2k=i \Fk\- A large portion of the memory requirement is due to the fact that, in order 
to achieve robustness, we need to compute the a-thresholds with long-integer arithmetic. The 
integers involved can get fairly long because the corresponding expressions are rather involved 
(see expressions (5-3) through (5-5) in section 5.3). However, as soon as their correct order is 
determined, the exact values of the a-thresholds, which are long-integer rationals, are no longer 
needed. For rendering purposes, standard floating-point accuracy is certainly sufficient. Duplicate 
a-thresholds, including the ones caused by attached simplices, are eliminated in the sorting phase. 
The size of the a-spectrum, that is, the final number of a-thresholds, plus the values and oo, is 
usually considerably smaller than the number of simplices in the triangulation. 



9 As a matter of fact, the second author implemented a variant of the randomized algorithm, achieving perfor- 
mance numbers roughly twice as fast than the ones in table 2 (see [34]). 
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Table 2: Performance of the incremental-flip algorithm. 













SoS primitives 


So 


S depth 


long-integer 


arithmetic 


memory 


CPU 


s 


\s\ 


\V\ 


|spectrum| 


Expr (5-5) 


Expr (5-4) 


Expr (5-3) 


max 


mean 


+ /" 


* 


Mb 


seconds 


molecule 


318 


8317 


4067 


1984 


1469 


632 





0.000000 


956974 


732611 


0.930 


7.80 


mal 


450 


11141 


5164 


2614 


1754 


1071 


37 


0.016009 


1232528 


968669 


1.172 


9.37 


tori 


800 


25193 


7811 


6003 


3338 


1848 





0.000000 


2808806 


2153544 


2.680 


27.23 


gridl 


1000 


25383 


5 


5961 


4854 


2700 


39 


0.490012 


2141341 


1911438 


2.562 


16.12 


spiral 


1248 


39849 


7939 


9556 


4710 


2649 


37 


0.009140 


4289522 


3357255 


4.195 


40.87 


molecule2 


1366 


37313 


20475 


8953 


8050 


3473 





0.000000 


4573791 


3592201 


4.141 


44.90 


universe 


1717 


44358 


21443 


10623 


7155 


3690 





0.000000 


5166676 


4077980 


4.539 


41.19 


pdb9pap 


1908 


52543 


31146 


12633 


12515 


5998 





0.000000 


5298887 


6665298 


5.672 


58.07 


ma2 


2206 


95056 


32434 


22779 


6611 


3046 





0.000000 


10081520 


7799049 


19.078 


82.66 


bust 


2630 


73021 


34041 


17475 


11799 


7541 


2 


0.006083 


8408574 


6776499 


8.023 


90.05 


f ract2 


2704 


73109 


31336 


17532 


9399 


7465 


24 


0.044352 


8010544 


6575623 


7.812 


83.62 


br 


3762 


106879 


61552 


25736 


24711 


11121 





0.000000 


13511423 


10827946 


11.383 


116.29 


liss_5_8 


4200 


121171 


56730 


28977 


17121 


11943 





0.000000 


14056712 


11096304 


12.781 


139.64 


teapot 


4668 


115725 


15434 


27142 


16626 


10366 





0.000000 


12494168 


10167622 


12.164 


128.27 


ku2 


5292 


154963 


60715 


37188 


33551 


19378 


39 


0.093015 


18000791 


15529705 


15.977 


153.95 


phone_l 


6070 


166331 


80590 


40014 


25444 


19773 


24 


0.000936 


19477178 


15850522 


18.234 


227.27 


rat_T2a 


9600 


261465 


36543 


62836 


56180 


28668 


24 


0.016854 


30276256 


24180721 


25.812 


237.65 


ra 


10000 


270343 


132353 


64804 


42307 


25245 


2 


0.000042 


32007436 


25724425 


28.445 


331.43 


hsr 


10088 


254241 


5763 


60974 


49469 


24431 


21 


0.003192 


23493214 


18959384 


24.836 


186.20 


smoke015 


15000 


400795 


200163 


96114 


65121 


38926 





0.000000 


47929411 


38659261 


43.898 


594.94 



Table 3: Performance of the a-interval generator. 



7 Applications and Further Illustrations 



It is important to point out that a-shapes are a fairly generic tool that can be used in many 
applications that have to do with shape, including automatic mesh generation and geometric 
modeling (see figure 3). Indeed, they can be used as a concrete expression of shape, which is 
often all that is needed. Similarly, three-dimensional a-shapes can be used to identify clusters 
in trivariate data. Beyond these generic applications, there are others that rely on particular 
properties of a-shapes. For these applications, it would be difficult to replace a-shapes by any 
other reasonable notion of shape. Two such applications are briefly addressed. 

Molecular Structures. Molecules are usually modeled as conglomerates of atoms with fixed 
relative positions. Each atom is represented by a ball around a center point, and the radius of the 
ball depends on what the model is supposed to express. For example, in the so-called space-filling 
diagram (see section 3.1) the balls encompass the idealized locations of the electrons so that balls 
of nearby atoms typically overlap. This diagram, defined as a union of balls, is in a strict geometric 
sense dual to the a-shape of the center points, assuming each ball has radius equal to a. The 
a-shape can thus be used to compute structural properties of the space-filling diagram, such as 
its connectivity. Alternatively, the a-shape itself, for this value of a, can be used to model and 
manipulate the molecule. When different atoms are represented by balls of different sizes then 
weighted a-shapes need to be used (see section 3.4). 

Molecules with intersting a-shapes arise in the study of proteins and how they fold (see figure 4). 
The geometric locations of the atoms of about 500 proteins are known today. However, there are 
many more gene sequences that can be determined. One of the goals of theoretical molecular biol- 
ogy is to obtain three-dimensional positional information from the knowledge of these sequences. 
This is the problem of protein folding [20,36]. Since the a-shape is computationally inexpensive 
and because it closely reflects the physical reality of molecules, it is hoped that a-shapes prove to 
be a useful tool in future protein folding research. 

Distribution of a Point Set. An interesting though ill-defined geometric problem arises in 
the study of the distribution of galaxies in our universe. As observed in studies like in [4,19], 
the galaxies are distributed in an unexpected and rather nonuniform manner (see figure 5). As- 
tronomers have measured the location of about 170000 galaxies, each one represented by a point 
in three-dimensional space. It appears that a large number of galaxies are located on or close 
to sheet-like and to filament-shaped structures. In other words, large subsets of the points are 
distributed in a predominantly two- or one-dimensional manner. 

How can this intuitive notion of the dimension of a point distribution be captured? A possible 
answer can be given by considering the entire family of a-shapes. Let A(a) be the surface area 
of the a-shape and let V(a) be its volume. To measure the degree to which the points are two- 
dimensionally distributed, it might be interesting to investigate the relative variation of A and V 
over the range of a values between and oo. More generally, it would be interesting to study 
the relationship between a-shapes ant the notion of fractals and fractal dimension. Through the 
availability of signatures of measures (see section 8), the a-shape gives comparably efficient access 
to metric information over a range of detail monitored by a. Resulting quantifications can be 
useful in the comparative study of the actually observed galaxy distribution and simulated data 
(see [10] for first steps in this direction). 
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Figure 3: Phone, [n = 6070, |F 2 | = 80131] 
The points are obtained from a public domain data set for modeling and rendering programs. All 
connectivity information (edges and triangles) is removed, and in order to generate more points, 
the centroids of all triangles are added. As described in section 3, each a-shape is triangulated by 
the tetrahedra of the corresponding a-complex. This might be useful in the automatic generation 
of meshes for objects with nonconvex surfaces. 
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Figure 4: Molecule, [n = 318, |F 2 | = 4000] 
The data represents a time-averaged molecular dynamics structure of gramicidin A, a peptide that 
forms a channel for ion and water movement across lipid membranes. The major structural motif 
is a right-handed beta-bonded helix. The tunnel of the macro-structure can be detected using 
relatively large a values. Smaller values of a result in a-shapes with larger numbers of isolated 
triangles and edges. These a-shapes reveal the helix of the micro-structure. 
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Figure 5: Universe, [n = 1717, |F 2 | = 21321] 
This data represents a simulation of the positions of galaxies within our universe. The theory is 
that galaxies first clustered into sheet-like structures, then progressed to filament-shaped structures 
at the intersection of multiple sheets. As filaments began to intersect, global clusters appeared. It 
is interesting to investigate the macro- and micro-structure of the galaxies, including the detection 
of large voids and local or global clusters. The full spectrum of a-shapes promises to be useful in 
this study. 
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8 Summary and Open Problems 



The main contribution of this paper is the introduction of a sound framework that formally 
captures the rather intuitive notion of the "shape" of a point set in space. This is the concept of 
three-dimensional a-shapes. A prototype version of a robust a-shape tool has been implemented. 
The authors of this paper hope that this tool will find many users within the engineering and 
the scientific computing and visualization communities. However, there is still a lot of work to 
be done. For example, the extensions mentioned in section 3.4 are worthwhile implementing, and 
this is planned in the near future. The extensions mentioned below are either less specific or 
theoretically not well understood. 

Improving the Running Time. A large fraction of the time used to construct a-shapes is needed 
for computing the Delaunay triangulation of the points. The algorithm used in our implementation 
takes time 0(n 2 ) in the worst case, independent of the number of simplices of T>. However, it 
rarely exhibits worst-case behavior. Still, it would be useful to have an algorithm whose running 
time is roughly proportional to the size of T>. Is it possible to construct T> in time O(n\og n + m), 
where m = \F\ U F 2 U F 3 \? A first step towards such an algorithm is the output-sensitive convex 
hull algorithm of Seidel [37]. If combined with the methods of Matousek and Schwarzkopf [30] it 
runs in time (9(n 4 / 3+e + mlogn) for n points in 1R 3 . By randomization, one can also achieve an 
expected running time roughly proportional to the expected number of simplices, if an underlying 
distribution is assumed. This is explained in [16]. 

On a different level, the running time of our program can be improved by speeding up the geometric 
primitives which all reduce to integer computations (see section 6). According to our experimental 
studies, about 75% of the time is spent doing integer arithmetic. This implies that appropriate 
hardware support might have a significant impact on the running time. 

Maintaining Alpha Shapes. In some applications it is necessary to construct a-shapes across 
a number of different point sets, and often these point sets are very similar to each other. For 
example, a point set might undergo local changes within an iterative process, and the a-shape or 
some feature of it is to be constructed at each step of the iteration. A local change might be the 
insertion of a new point, the deletion of an old point, the dislocation of one point, the dislocation 
of a subset of the points, etc. The development of efficient update algorithms that reuse available 
structure as much as possible can lead to dramatic improvements of the overall running time. 

Features and Signatures. Individual a-shapes are interesting geometric objects, and it would 
be useful to have efficient algorithms that can analyze its geometric and topological properties 
or features. For example, computing the volume is fairly straightforward because the a-complex 
provides a triangulation of the a-shape, and the volume of the a-shape is simply the sum of the 
volumes of the tetrahedra. More challenging is the computation of the connectivity of the a-shape 
as expressed by its first three homology groups (see for example Giblin [21]). 

As suggested by the second application discussed in section 7, the history of a feature, over all 
values of a from through oo, is of interest. Consider some specific feature, say, the number of 
connected components of the a-shape. The corresponding signature is a function c: [0, +oo] — > M } 
so that c(a) is the number of components of S a . This function reflects the evolution of the number 
of components as a changes continuously from to +oo. Given the a-spectrum, it is fairly easy to 
compute c. Start at a = and maintain a union-find data structure (see for example [5]) storing 
the components as threshold values are processed in increasing order. A more challenging task is 
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the computation of the signatures for higher-order homology groups. Such signatures might be 
handy in the selection of an appropriate a value which typically depends on the application and 
the user's momentary focus of attention. 



Available Software 

The implementation mentioned in section 6 is called "Alvis — A 3D Alpha Shape Visualizer." 
This tool allows the user to interactively select a-values and display the corresponding shape. 
A small collection of signatures aids the selection process. We see the program as an extension 
to the paper, as an "animated figure," so to speak. One of its purposes is to effectively convey 
the concept of a-shapes to the engeneering and scientific computing community. It is available 
via anonymous ftp from ftp.ncsa.uiuc.edu (141.142.20.50). The latest release of the code, 
including the code for three-dimensional Delaunay triangulations, can be found in the directory 
SGI/Alpha-shape. The code requires a Silicon Graphics workstation, running Irix 4.0, or later; 32 
Mb main memory, or more, are advisable. This code is a new kind of shareware: you share with 
us your experience in applying Alvis to engineering and science problems, and we do our best to 
develop the software so it can meet your needs. For questions contact <alpha@ncsa .uiuc . edu>. 
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"Don't look like a convex hull . . . get yourself in a-shape!" 

— David Knapp. 
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